 ;;########################################################################
;; missd7.lsp
;; Copyright (c) 1998 by Pedro Valero (valerop@uv.es)
;; Code for multiple imputation of Missing data 
;; 
;; 
;;########################################################################

(defmeth missing-data-model-object-proto :da-options 
  (&key (dialog t) (num-iterations 100) (num-datasets 3) (sample-every-k-iterations 10))
  (when (send self :dialog)
        (let* ((result nil)
               (iterations-text (send text-item-proto :new 
                                      "Number of Iterations (minimum=110)"))
               (num-datasets-text (send text-item-proto :new 
                                        "Number of datasets (you can add later)"))
               (sample-iterations-text (send text-item-proto :new 
                                        "Number of iterations between files"))
               (iterations (send edit-text-item-proto :new " 1000"))
               (datasets (send edit-text-item-proto :new "   3"))
               (sample-iterations (send edit-text-item-proto :new "  30"))
               (num-iterations)
               (num-datasets)
               (sample-every-k-iterations)
               (Ok (send modal-button-proto :new "OK" 
                         :action #'(lambda ()
                                     (list
                                      (setf num-iterations
                                            (read-from-string 
                                             (send  iterations :text)))
                                      (setf num-iterations
                                            (if (< num-iterations 110) 110 num-iterations))
                                      (setf num-datasets
                                            (read-from-string (send datasets :text)))
                                      (setf sample-every-k-iterations
                                            (read-from-string (send sample-iterations :text)))
                                      ))))
               (cancel (send modal-button-proto :new "Cancel" 
                             :action #'(lambda () 
                                         (send dialog :modal-dialog-return nil))))

               (dialog    (send modal-dialog-proto :new
                                 (list (list iterations-text iterations) 
                                       (list num-datasets-text datasets)
                                       (list sample-iterations-text sample-iterations)
                                       (list ok cancel))))
               (a (send dialog :default-button ok))
               (input-options (send dialog :modal-dialog))
               )
  (send self :start-da num-iterations num-datasets sample-every-k-iterations))))

(defmeth missing-data-model-object-proto :start-da (num-iterations num-datasets sample-every-last-k)
  (let* (
         (datasets num-datasets)
         (emcovariance (send self :emcovariance))
         (em-means (combine (send self :em-means)))
         (iter num-iterations)
         (data (send self :data))
         (res nil)
         (output-dataset (mapcar #'(lambda (val) 
                                     (- iter (* val sample-every-last-k)))
                                 (iseq 1 num-datasets)))
         )
    (send self :data-augmentation  
          iter
          emcovariance
          em-means
          data
          output-dataset)))


  (defmeth missing-data-model-object-proto :add-dataset-multiple-imputed-list (dataset)
    (send self :multiple-imputed-datasets-list
          (append (list dataset) 
                  (send self :multiple-imputed-datasets-list))))


  

  (defmeth missing-data-model-object-proto :data-augmentation (iter covariances means data output-dataset)
    (let* ( 
            (prev-means (mapcar #'(lambda (col) 
                                    (mean (non-missing col))) 
                                (column-list data)))
            (prev-st (mapcar #'(lambda (col) 
                                 (standard-deviation (non-missing col))) 
                             (column-list data)))
            (mat-st (make-array (list 1 (length prev-st)) :initial-contents prev-st))
            (mat-st (matmult (transpose mat-st) mat-st))
            (orig-covar)
            (orig-means);these variables are related with the series of the procedure
            (orig-matrix)
            (data-orig (copy-array data))
            (data (copy-array (normalize (center data))))           
            (res-I-step nil)
            (res-p-step nil)
            (iter iter)
            (dim-arrays (length means))
            (covariances (identity-matrix dim-arrays))
            (means (repeat '0 dim-arrays))
            (series (make-array (list (1+ dim-arrays) (1+ dim-arrays) iter)))
            (T-obs (first (tobs (copy-array data))))
            (res-p-step (list covariances means))
            (n (array-dimension data 0))          
            (m-0 -1)
            (V-0 (make-array (list dim-arrays dim-arrays) 
                             :initial-element '0))
            (r-0 0)
            (u-0 (repeat 0 dim-arrays))
            (m-p nil)
            (v-p nil)
            (u-p nil)
            (r-p nil)
            (patterns  (patterns-missing data))
            (rows-in-patterns (first (cases-in-missing-patterns (copy-array data))))
            (hyper-parametters nil)
            (missing-in-missing-pattern-list (Missing-in-missing-pattern-list patterns))
            (observed-in-missing-pattern-list (observed-in-missing-pattern-list patterns))    
            (sw-res-i-step nil)
            (orig-data nil)
            (output-dataset output-dataset)
            (creator (send *desktop* :selected-icon))
            (start-time)
            (time-iteration)
            )
      (setf *multiple-imputation-report-window* 
            (report-header "Missing Data Multiple Imputation"
                           :scroll t :location '(250 100)))
      ;(send *multiple-imputation-report-window* :top-most?) This is not the message
      (add-text *multiple-imputation-report-window* 
                (format nil "Multiple Data Imputation") :scroll t)
      (add-text *multiple-imputation-report-window* 
                (format nil "~%Copyright (c) 1998-2000 Pedro Valero (valerop@uv.es)") :scroll t)
      (add-text *multiple-imputation-report-window* 
                (format nil "~%Multiple Imputation can be a lengthy process") :scroll t)
      (add-text *multiple-imputation-report-window* 
                (format nil (strcat "~%You have selected " (princ-to-string iter) " iterations")) :scroll t)
      (add-text *multiple-imputation-report-window* 
                (format nil "~%Notice that after the multiple imputation process finishes") :scroll t)
       (add-text *multiple-imputation-report-window* 
                (format nil "~%a new visualization is available in the missing data object.") :scroll t)
      
      (time (dotimes (i iter)
                     (cond 
                     
                     ;this part prints an estimation of the time for the computation
                            ((= i 10)
                             (setf start-time (get-internal-real-time)))
                            ((= i 30)
                             (setf time-iteration (/ (/ (- (get-internal-real-time)
                                                           start-time)
                                                        internal-time-units-per-second) 20))
                             (add-text *multiple-imputation-report-window* 
                                       (format nil "~%Time to finish aprox. (in seconds): " :scroll t))
                             (add-text *multiple-imputation-report-window*  
                                       (format nil (princ-to-string (round (* iter time-iteration))))))
                            ((= i (- iter 1))
                             (add-text *multiple-imputation-report-window* 
                                       (format nil "~%Done" :scroll t)
                                       (send *multiple-imputation-report-window* :show-window)))
                           
                            ((and (> i 30)
                                (= (denominator (rational (/ i 50))) 1))
                             (add-text
                              *multiple-imputation-report-window*  
                              (format nil 
                                      (princ-to-string (strcat "~%The time left is:" 
                                                               (princ-to-string (round (* (- iter i) time-iteration)))
                                                               "(s) after "(princ-to-string i) " iterations")))))
                            )
                                                             
                           
                     
                     ;end of time computation
                           
               (setf res-i-step 
                     (send self :I-step 
                           (first res-p-step)    
                           (second res-p-step)
                           (copy-array data)
                           (copy-array t-obs)
                           patterns
                           rows-in-patterns
                           missing-in-missing-pattern-list
                           observed-in-missing-pattern-list
                           ))
               (setf sw-res-i-step (second res-i-step))
               (setf (aref sw-res-i-step 0 0) 1)
               (setf sw-res-i-step 
                     (first (schafer-sweep-operator 
                             sw-res-i-step (list 0))))
               


               #| (setf means-orig (+ prev-means
                                      (* prev-st (second res-p-step))))
               (push means-orig b)
               (push (second res-p-step) c)|#


                     (setf means (combine (select sw-res-i-step 0 
                                            (iseq 1 (1- 
                                                     (array-dimension 
                                                      (second res-i-step) 0))))))
               
                     (setf covariances (select sw-res-i-step
                                         (iseq 1  (1- (array-dimension 
                                                       (second res-i-step) 0)))
                                         (iseq 1  (1- (array-dimension 
                                                       (second res-i-step) 0)))))
                     (setf orig-covar (matmult 
                                       mat-st
                                       covariances))
                     (setf orig-means (+ (* prev-st means) prev-means))
                     (setf orig-matrix (my-border-matrix orig-covar orig-means orig-means n))
                           
                     (setf (select series 
                                   (iseq (1+ dim-arrays))  
                                   (iseq (1+ dim-arrays)) i)
                           (make-array (list (1+ dim-arrays) (1+ dim-arrays) 1) 
                                       :initial-contents (combine orig-matrix)))
               
               (setf hyper-parametters
                     (send self :update-hyperparametters 
                           n 
                           m-0
                           v-0
                           u-0
                           r-0
                           means
                           covariances))
               (setf m-p (first hyper-parametters))
               (setf V-p (second hyper-parametters))             
               (setf u-p (third hyper-parametters))
               (setf r-p (fourth hyper-parametters))
               (setf res-p-step                 
                     (send self :P-step  n m-p v-p u-p r-p))
               (when (member i output-dataset) (send self :create-files res-i-step prev-means prev-st creator))
                    (send self :da-iterations-info series)
               ))
      (list means covariances orig-data)
      ))
  
   (defmeth missing-data-model-object-proto :create-files (res-i-step prev-means prev-st creator)
     (let* (
            (res-i-step res-i-step)
            (prev-means prev-means)
            (prev-st prev-st)
            (creator creator)
            (orig-data
             (apply 'bind-columns 
                    (mapcar #'(lambda (col-orig col m st)
                                (let ((miss-in-var 
                                       (which (mapcar #'(lambda (val) (equal val 'nil))
                                                      (combine col-orig)))))
                                  (when miss-in-var (setf (select col-orig miss-in-var)
                                                          (+ (* (select  col miss-in-var) st) m)))
                                  col-orig))
                            (column-list (send self :data))
                            (column-list (third res-i-step)) 
                            prev-means prev-st))))
       (send self :add-dataset-multiple-imputed-list
             (data (send self :name)     
                   :created  creator
                   :title     (concatenate 'string "MI-"(send self :title))
                   :variables (send self :variables)  
                   :data     (combine orig-data)
                   :labels (send self :labels)
                   ))))


  (defmeth missing-data-model-object-proto :update-hyperparametters 
  (n m-p V-p u-p r-p means covariances)
    (let* (       
                  (V-p 
                   (+  v-p
                       (* n covariances) 
                       (matmult (/ (* r-p n) (+ r-p n)) 
                                (matmult (make-array 
                                          (list  (length means) 1 ) 
                                          :initial-contents (combine (- means u-p  ))) 
                                         (transpose (make-array (list  (length means) 1 ) 
                                                                :initial-contents (combine 
                                                                                   (- means u-p))))))))
                  (u-p (+ (* (/ n (+ r-p n)) means) (* (/ r-p (+ r-p n)) u-p)))
                  (r-p (+ r-p n))
                  
                  (m-p (+ m-p n)))       
      
      (list m-p  V-p u-p r-p)
      ))
  
  (defmeth missing-data-model-object-proto
    :I-step 
    (covariances means data t-obs patterns rows-in-patterns missing-in-missing-pattern-list observed-in-missing-pattern-list) 
    (let* 
      (
       (patterns patterns)
       (rows-in-patterns rows-in-patterns)
       (y-ij (copy-array data))
       (parametters (my-border-matrix 
                     covariances (combine means) (combine means)
                     -1))    
       (c (copy-array parametters))
       (s (length rows-in-patterns))
       (p (- (array-dimension parametters 0) 1))    
       ; (c (make-array (list p 1)
                        ;               :initial-element 0))
       (n-cases (array-dimension data 0))
       (r (make-array (list s p) :initial-contents patterns))
       (t-obs (copy-array t-obs))
       (C-miss (make-array (array-dimensions parametters)))
       (z-simulated nil) 
       (product-tobs nil)
       (suma-tobs nil)
       )
      (dotimes (i s)       
               (dotimes (j p)
                        (if (and (equalp (aref r i j) 1) 
                                 (> (aref parametters (+ 1 j) (+ 1 j)) 0))
                            
                            (setf parametters 
                                  (select 
                                   (schafer-sweep-operator 
                                    parametters (list (+ 1 j))) 0))
                            )                           
                        (if  (and (equalp (aref r i j) 0) 
                                  (< (aref parametters (+ 1 j) (+ 1 j)) 0))                              
                             (setf parametters 
                                   (select  
                                    (reverse-schafer-sweep-operator 
                                     parametters (list (+ 1 j))) 
                                    0))))
               
               (when (select missing-in-missing-pattern-list i)
                     (setf (select C-miss 
                                   (1+ (select missing-in-missing-pattern-list i))
                                   (1+ (select missing-in-missing-pattern-list i)))
                           (first (chol-decomp
                                   (select parametters
                                           (1+ (select missing-in-missing-pattern-list i))
                                       (1+ (select missing-in-missing-pattern-list i))
                                           ))))
                     
                      #|(dolist (m (select rows-in-patterns i))
                             (dolist (n (select missing-in-missing-pattern-list i))     
                                     (setf (aref y-ij m n) (aref parametters   (1+ n) 0))  
                                     (dolist (k (select observed-in-missing-pattern-list i))
                                             (setf (aref y-ij m n) 
                                                   (+ (aref y-ij m n)
                                                      (* (aref parametters (1+ k) (1+ n))
                                                         (aref y-ij m k)))))
                                     
                                 (dolist (k (select missing-in-missing-pattern-list i))
                                         (when (<= k n)
                                               (setf (aref y-ij m n) 
                                                     (+ (aref y-ij m n)
                                                        (* (aref c-miss (1+ k) (1+ n))
                                                           (first (normal-rand 1)))
                                                        ))))
                                   ))|#
                     
                     (setf (select y-ij (select rows-in-patterns i) 
                                   (select missing-in-missing-pattern-list i))
                           (apply 'bind-columns
                                  (mapcar #'(lambda (val col) (+ val col))
                                          (combine (select parametters (1+ (select missing-in-missing-pattern-list i)) 0))
                                          (column-list (transpose
                                                        (matmult 
                                                         (select parametters (1+ (select missing-in-missing-pattern-list i))
                                                                 (1+ (select observed-in-missing-pattern-list i)))
                                                         (transpose (select y-ij (select rows-in-patterns i) 
                                                                            (select observed-in-missing-pattern-list i)))))))))
                     (setf (select y-ij (select rows-in-patterns i) 
                                   (select missing-in-missing-pattern-list i))
                           (apply 'bind-columns
                                  (map-elements #'+
                                                (*  (diagonal (select c-miss 
                                                                      (1+ (select missing-in-missing-pattern-list i))
                                                                      (1+ (select missing-in-missing-pattern-list i))))
                                                    (normal-rand
                                                     (length (select missing-in-missing-pattern-list i))))
                                                (column-list (select y-ij 
                                                                     (select rows-in-patterns i) 
                                                                     (select missing-in-missing-pattern-list i))))))

                             ))
                              
                              (setf suma-tobs (map-elements 'sum (column-list y-ij)))
                              (setf product-tobs (cross-product y-ij))
                              (setf t-obs (my-border-matrix product-tobs suma-tobs suma-tobs (array-dimension y-ij 0)))
                              
                              (list t-obs  
                                    (first (schafer-sweep-operator 
                                            (* (/ 1  (array-dimension y-ij 0))
                                               (copy-array t-obs)) 
                                            (list 0)))
                                    y-ij )
                              ))
      
      
      
(defmeth missing-data-model-object-proto :p-step  (n m V u-0 r)
  (let* (
         (m m)
         (p (array-dimension V 0))
               (v v)
         (chol-V (first (chol-decomp V)))         
         (b-ji (make-array (list p p) :initial-element 0))
         (res-M nil)
         (u-0 u-0)
         (r r)
         )    
    
   (dotimes (i p)
             (dotimes (j p)
                      (when (> j i)
                            (setf (select b-ji i j) 
                                  (first (normal-rand 1))))
                      (when (= j i)
                            (setf (select b-ji j i) 
                                  (first (sqrt (chisq-rand 1 (+ 1 (- m (1+ j))))))))))
    
    (setf res-M (solve (transpose b-ji)  (transpose chol-v)))
    (list (cross-product res-M) 
          (+ u-0 (* (/ 1 (sqrt r)) 
                    (matmult   (transpose res-m) (normal-rand p)))))
    ))      
 

#|(defmeth missing-data-model-object-proto :p-step  (n m V u-0 r)
  (let* (
         (m m)
         (p (array-dimension V 0))
               (v v)
         (chol-V (first (chol-decomp V)))         
         (b-ji (make-array (list p p) :initial-contents (normal-rand (* p p))))
         (res-M nil)
         (u-0 u-0)
         (r r)
         )    
    (mapcar #'(lambda (j) 
                (setf (select b-ji j j)  (first(sqrt (chisq-rand 1 (+ 1 (- m (1+ j))))))))
            (iseq p))
    (setf res-M (solve (transpose b-ji)  (transpose chol-v)))
    (list (cross-product res-M) 
          (+ u-0 (* (/ 1 (sqrt r)) 
                    (matmult   (transpose res-m) (normal-rand p)))))
    ))|#
          
    


(defmeth missing-data-model-object-proto :visualize-da-imputation ()
  (let* (
         (ct (make-container :free t :type 2 :local-menus t :show nil))
         (temp   (enable-container ct))
         (list-par (combine (mapcar '(lambda (x)
                                       (concatenate 'string "mean_" x ))
                            (send self :variables))
                            (mapcar '(lambda (x)
                                       (concatenate 'string "Covar_" x ))
                                    (first (send self :listvariables)))
                            (mapcar '(lambda (x)
                                       (concatenate 'string "Variance_" x ))
                                    (send self :variables))))
         (aref-means (mapcar #'(lambda (val) (combine 0 val))
                             (iseq 1 (length (send self :variables)))))
         (aref-variances (mapcar #'(lambda (val) (combine val val))
                                (iseq 1 (length (send self :variables)))))
         (aref-covariances (mapcar #'(lambda (val1 val2) (combine val1 val2))
                                   (1+ (second (send self :listvariables))) 
                                   (1+ (third (send self :listvariables)))))
         (arefs (append aref-means aref-covariances aref-variances))
         (list-vars 
          (name-list  
           list-par :show nil :title "Parameters"))
         (num-iter (array-dimension (send self :da-iterations-info) 2))
         (h 
          (combine
           (select (send self :da-iterations-info) 0 1
                   (iseq (third (array-dimensions (send self :da-iterations-info)))))))
         (ac (autocovariances h))
         (plines (plot-lines (iseq (length h)) h :color 'blue :title (strcat "Time Series for " (select list-par 0)) ))
         (ac-plot (plot-lines (iseq (length ac)) ac :color 'blue :title (strcat "Autocorrelation for " (select list-par 0))))
         (missing-model-object self)
         )
  

    (spreadplot (matrix (list  2 3) 
                        (list list-vars ac-plot nil nil plines nil))
                :container ct 
                :span-right (matrix (list 2 3) (list 1 2 0 1 2 0))
                :span-down (matrix (list 2 3) (list 2 1 0 0 1 0)) :show t)
    (send ac-plot :legend1 "Data Augmentation series")
    (send plines :legend1 "Data Augmentation series")
    (send plines :redraw)
    (send ac-plot :range '1 -1 1)
    (send ac-plot :add-lines (send ac-plot :range '0) (list 0 0) :color 'orange :type 'dashed)
    (send ac-plot :add-lines (iseq (length ac)) ac :color 'blue)
    (defmeth list-vars :do-click (a b c d)
      (call-next-method a b c d)
      (setf selection (first (combine (send self :selection))))
      (setf h (combine (select (send missing-model-object :da-iterations-info) 
                          (mapcar #'(lambda (val) (first val)) (list (select arefs selection)))
                          (mapcar #'(lambda (val) (second val)) (list (select arefs selection)))
                          (iseq num-iter))))
      (send plines :start-buffering)
      (send plines :clear-lines)
      (send plines :clear-points)
      (send plines :legend2 (strcat "Time Series for " (select list-par selection)))

      (send plines :add-lines (iseq num-iter) h :color 'blue)
     ; (send plines :add-points (iseq num-iter) h :color 'blue :symbol 'dot)
            
      (send plines :adjust-to-data)
      (send plines :buffer-to-screen)
      (send ac-plot :start-buffering)
      (setf ac (autocovariances h))
      (send ac-plot :clear-lines)
      (send ac-plot :clear-points)
      (send ac-plot :legend2 (strcat "Autocorrelation for "  (select list-par selection)))
      (send ac-plot :add-lines (send ac-plot :range '0) (list 0 0) :color 'orange :type 'dashed)
      (send ac-plot :add-lines (iseq (length ac)) ac :color 'blue)
      
     ; (send ac-plot :add-points (iseq (length ac)) ac :color 'blue :symbol 'cross)
      (send ac-plot :redraw)
      (send ac-plot :buffer-to-screen)

      )
    (disable-container)
    ))

(defun autocovariance (data lag &key (type "biased"))
    (let* ((n (length data))
           (m (mean data))
           (z (- data m))
           (c (abs lag))
           (d (if (string= type "unbiased") (- n c) n))
           (a (if (string= type "circular") z (butlast z c)))
           (b (if (string= type "circular") (shift-sequence z c) (butfirst z c))))
      (if (<= c n) (/ (inner-product a b) d) 0.0)
      ))

(defun autocovariances (var &optional (numlags 100))
  (if (variancep var)
     (let*
      (
       (var (/ (- var (mean var)) (standard-deviation var)))
       (autocovariances (mapcar #'(lambda (lag) (autocovariance var lag :type "unbiased")) (iseq 1 numlags)))
       )
       autocovariances)
      (repeat '0 numlags)))
  
  (defun autocovariance-plot (var &optional (numlags 100))
    (let*
      (
       (var (/ (- var (mean var)) (standard-deviation var)))
       (autocovariances (mapcar #'(lambda (lag) (autocovariance var lag :type "unbiased")) (iseq 1 numlags)))
       )
      (plot-lines (iseq 1 numlags) autocovariances)
      ))
  
  (defun butfirst (x &optional (n 1))
    (select x (which (<= n (iseq (length x))))))


  #|(defun test ()
  (setf b nil)
  (setf c nil)
  (send current-model :mi)
  
  (setf h (reverse (first (last (transpose b)))))
  (setf c (select (transpose c) 10))
  (setf 

  (spreadplot (matrix (list  1 2) 
                      (list (autocovariance-plot h) 
                            (plot-lines (iseq (length h)) h))) :container ct
              :show t)
  (disable-container)
  
  )
|#